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ABSTRACT 

We study the influence of gas metaflicity, turbulence, and non-equilibrium chemistry on the 
evolution of the two-phase interstellar medium (warm and cold atomic phases), and thereby 
constrain the initial conditions for star formation prevailing in turbulent gas. We perform high- 
resolution simulations in three dimensions, including a realistic non-equilibrium treatment of the 
ionization state of the gas, and examine both driven and decaying turbulence. This allows us 
to explore variations in the metallicity Z. In this paper, we study solar metallicity, Z = Zq, 
and low metallicity, Z = 10~^ Zq, gas. For driven, large-scale turbulence, we find that the 
influence of the metallicity on the amount of mass in the cold gas component is small. However, 
in decaying turbulent conditions this picture is much changed. While cold regions survive in the 
case of solar metallicity, they are quickly heated and dispersed in low-metallicity gas. This result 
suggests that star formation can be suppressed in environments of low metallicity, unless a strong 
turbulent driver is acting on time scales shorter than a few turbulent crossing times. Inter alia 
this finding could explain the overall ineflicient star formation as well as the burst-like mode of 
star formation found in metal-poor, gas-rich systems like dwarf galaxies. 

Subject headings: galaxies: high-redshift - galaxies: ISM - galajcies: star formation - ISM: clouds - ISM: 
evolution - ISM: general 
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1. Introduction 

The cold interstellar medium (ISM) and, in par- 
ticular, the molecular part of it, forms the reser- 



voir 


available for star formation (Wong & Blitz 


2002 




Kennicutt et al. 


2007 




Bigiel et al. 


20081. 



The properties of this cold gas component regulate 
the mass distribution of prestellar cores in a molec- 



ular cloud (e.g., [Padoan fc Nordlund||2002 |2011| 
Hennebelle fc Cha brier'2008', 'Heits ch et al.[2008j , 
and hence the stellar initial mass function (IMF; 
e.g.,|Motte et al.|1998| [Nutter fc Ward-Thompson] 
2007). It has been suggested that star formation 



within a molecular cloud is regulated by the inter- 



play of gravity and turbulence (e.g., Klessen et al. 
2000 1 . However, on galactic scales, understand- 



ing the formation and the properties of the cold 
gas reservoir available for star formation is equally 
important. In addition to turbulence and grav- 
itational collapse, heating and cooling processes 
essentially determine the distribution of the cold 
ISM and therefore the global star formation effi- 



ciency ( [Ostriker et al.j|2010J . 

The ISM is highly turbulent, with gas veloc- 
ity dispersions of up to 50 times the sound speed 
of the molecular (10 K) gas component in local 
galaxies, and even higher factors in high-redshift 
disk galaxies (Genzel et al. 2006). However, tur- 



bulence has been shown to dissipate quickly ( Mac 



Low 1999) if the external driving is stopped, re- 



sulting in the need to continuously drive turbu- 



lence either locally via stellar feedback (Ostriker 
et al. 2010) or globally via, e.g., galaxy mergers 



( [Teyssier et"aI]|2010[|Karl et al.|2010|), clumpy in 



fall of cold gas ( [Elmegreen fc Burkert|2010 Aumer 
et al.|[20T0 ), or spiral density waves 

At the same time, the ISM is regulated by 
metallicity-dependent heating and cooling pro- 
cesses that tend to segregate it into distinct, ther- 
mally stable phases (warm and cold), which can 
coexist in approximate pressure equilibrium. The 
transition between the phases is driven by ther- 



mal instability (TI; |Field||1965| [Burkert fc Lin 
2000[), mo re specifically by its isobaric mode (e.g., 



PikerNer|[T968l |McKee fc Ostriker|[T977l |Wolfire| 



et al. 1995). Thus, for several reasons the effi- 
ciency with which a galaxy is able to form cold gas 
depends on the metal content of the gas. First, 
the turbulent shock compression of the ISM de- 
pends on the metallicity-dependent cooling rate 



(e.g., Omukai et al.|2005 ). Second, within the cold 
component, the dust abundance (i.e., metallicity) 
controls the conditions under which the transition 



from cold atomic to molecular gas occurs ( Gnedin 
et al.||2009l [Krumholz et al.|[2009l [Glover fc Mac 



Low 



20111. 



Third, depending on the strength 
of the interstellar radiation field, low gas metal- 
licities (< O.OIZq) alter the balance of heating 
and cooling in such a way that the development 
of a bi-stable phase can be suppressed (see semi- 
analytical estimates by [Spaans fc CaroUo 1997 



Spaans fc Norman|1997[ ) . This conclusion has also 
been drawn by [Schaye[ ([20041) who used equilib- 



rium models and the CLOUDY software to demon- 
strate that the formation of cold gas is a metallic- 
ity dependent process. His models are successfully 
employed to predict the outer edges of standard 
star forming galactic disks, for both regular spi- 
rals and dwarfs ( Leroy et al.[[2008 ). 



The idea that the (extra-)galactic star forma- 
tion efficiencies are connected to the present ISM 
properties has been further discussed in several 



recently developed theoretical models |Elmegreen 


20021 IKrumholz fc McKee,j2005| 


Schaye fc Dalla 


Vccchia 2008| [Gnedin et al.|2009 


Krumholz et al. 


2009f |Pelupessy fc Papadopoulos 


2009 Krumholz 


et al. 


2011). Some of them estimate the effects 



of turbulence and have shown that the thermody- 
namics of the ISM can control the resulting shape 
of the density probability density function (PDF) 



of the ISM gas ( 


Wada fc Norman[[2007 


Tasker fc 


Bryan 


2008 


Robertson fc Kravtsov 


20081. The 



numerical simulations we present here indicate 
that the density PDF is only strongly affected in 
the non-driven turbulent case. 

Previous highly resolved studies of the statis- 
tics of driven or decaying turbulence mostly as- 
sumed the gas to be isothermal (e.g.. 



Kritsuk 



eraLl[2007| [Beresnyak fc Lazarian[ [20091 includ 
ing magnetic fields). First steps toward a detailed 
numerical model of turbulence including gas ther- 
modynamics have been published by j^Kritsuk 



Norman (20041, Joung fc Mac Low (2006), Bran 



denburg et al. (2007), and de Avillez & Breitschw- 



|erdt^ (2007 ), who include parameterized cooling 
functions at solar metallicity. Kritsuk fc Nor-[ 



man 



(2004) and Brandenburg et al. (2007) study 



the amount of turbulence generated by the TI it- 



self, whereas Joung & Mac Low (2006) and de 



Avillez fc Breitschwerdt ( 2007 ) use supernova ex- 
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plosions as internal drivers of turbulence. Other 
authors have studied molecular cloud formation in 



colliding flows in two dimensions 


(Audit & Hen- 


nebelle||2005 


Heitsch et al.||2006 


2008 ) and three 


dimensions ( 


Slyz et al. 2005 Vazquez-Semadeni 


et al.||2007||Hennebelle et al.|2008 


Banerjee et al. 



2009) using parameterized cooling and heating 
processes. In these simulations turbulence is gen- 
erated via various instabilities occurring in the col- 
lision layer of the flow. These studies confirm the 
development of a multi-peaked density and tem- 
perature PDF when turbulence and equilibrium- 
based gas thermodynamics are treated simulta- 
neously. On the other hand, |Vazquez-SemadeiiI| 
(2009) conclude that turbulence acts faster than 



TI, and thus relegates classic instabilities (Field 



1965) to a secondary role in controlling the gas 
dynamics. He shows that shock-compressed dif- 
fuse gas can stay in an unstable regime for a few 
megayears before cooling efficiently. This is sup- 
ported by recent models of |de Avillez fc Breitschw- 



erdt (2010) who point out that the treatment of 
non-equilibrium ionization is essential for a correct 
treatment of the highly dynamical ISM evolution, 
at least in the high temperature regime. 

In this paper, we show that the formation and 
evolution of the cold ISM phase is a highly metal- 
licity dependent process. We base this conclusion 
on high-resolution, three-dimensional simulations 



with the Eulerian grid code FLASH (Fryxell et al. 



2000), including both driven and decaying turbu- 
lence, as well as metallicity-dependent heating and 
cooling ( Glover et al.||2010 ). We demonstrate the 
major influence of the gas metallicity in decaying 
conditions - a finding which highlights the global 
metallicity dependence of star formation. We sug- 
gest that this result may explain puzzling obser- 
vations, like the low star formation rates in metal- 
poor, gas-rich galaxies. 

2. Numerical Method 

We use FLASH 3.2, an MPI parallel Eulerian 
grid code, which solves the Euler equations for 
fluid flow. We modify the energy equation to in- 



clude heating and cooling effects (see Section 2.2 1. 
The implemented directionally split, finite differ- 
ence scheme is based on the piecewise parabolic 



method (PPM; [Woodward & Colella|1984 


Colella 


k Woodward 1984 


) in its direct Eulerian form. 



PPM allows one to model strong shocks and dis- 
continuities with high accuracy. Using the method 
of characteristics, the nonlinear flux of quantities 
between zones is obtained by solving a Riemann 
problem at each zone boundary in alternating one- 
dimcnsional sweeps through the grid. We use the 
hybrid Riemann solver (HLLE) inside shocks to 
avoid odd-even instabilities ( Hussaini et al.|1997 |. 



2.1. Turbulence 

At every time step, supersonic turbulence is 
driven on the largest spatial scales (1 < fc < 3, 
where fc = 1 corresponds to 27r/500pc and fc = 3 
corresponds to 27r/166pc). The driving energy 
input rate per unit mass, dEi^/dt = 1.04 x 
10~^erg g~^s~^, is chosen such that supersonic 
turbulence develops, with a roughly constant root- 
mean-squared (rms) velocity of 50kms~^. We use 
a mixture of solenoidal and compressive forcing 
with random phases and time-correlated driving 
forces, which leads to coherent large-scale struc- 
tures. The phases of the driving modes are mod- 
ulated on a turnover timescale, r^yn = 50 Myr, 
which is equal to five turbulent crossing times, or 
about 70% of a sound crossing time for gas at the 
initial temperature of 5000 K. For each mode and 
at each time step six separate phases (real and 
imaginary in three dimensions) are evolved by an 
Ornstein-Uhlenbeck random process, which repre- 
sents a velocity in Fourier space. 

2.2. Chemistry and Cooling 

Gas dynamics and gas chemistry are strongly 
coupled and should be solved at the same time 
( Glover et al.||20f0 ). In order to study the metal- 
licity dependence of the multi-phase ISM, we in- 
clude a basic and fast chemical network, based 



on Glover & Mac Low (2007), to follow a limited 
amount of non-equilibrium gas chemistry, specifi- 
cally the chemistry involved in the determination 
of the ionization fraction of the gas. When solving 
the energy equation we take account of radiative 
cooling due to line emission (Lya, Cn, Oi, and 
Sin) and dust, and the chemical cooling produced 
by the coUisional ionization and radiative recom- 
bination of hydrogen. We account for heating due 
to cosmic rays. X-rays, and photoelectric emission 
from small grains and polycyclic aromatic hydro- 
carbons (PAHs). For the ultraviolet background 
radiation field, we use the standard estimate of 
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Draine ( 1978 ). In this study, we neglect the effects 



of dust shielding. This is a reasonable approxima- 
tion in the warm gas in our solar metalhcity simu- 
lation, but may cause us to overestimate the tem- 
perature of the dense component; a colder dense 
component would in any case strengthen our con- 
clusions. It is a very good approximation through- 
out our low-metallicity simulation, as in this case 
the dust extinction is very small even in the cold 
clumps. For the cosmic ray ionization rate of neu- 
tral hydrogen, we adopt a value Ch = 3x 10~-^''s~^, 
while for the X-ray ionization and heating rates, 
we use the values given in Appendix A of |Wolfire| 
et al. (1995), assuming a fixed absorbing column 
of hydrogen Nh = 10^° cm~^. 

We do not include the formation and destruc- 
tion of molecular hydrogen (H2) in our current 
models. This may be a significant coolant in 



the warm, diffuse ISM in the Galaxy (Gnedin & 



Kravtsov 2010), although this remains unclear. 



However, it is very unlikely to be a significant 
coolant in the warm gas in our low-metallicity 
simulation, as its equilibrium abundance would 
be very small, even if gas-phase formation mech- 
anisms were included (Glover 2003). We would 



therefore expect the inclusion of H2 chemistry 
and cooling to merely heighten the differences be- 
tween the results of our solar metalhcity and low- 
metallicity runs. 

The heating and cooling terms, and also the 
chemical source and sink terms, are operator split 
from the advection and pdV work terms and are 



solved implicitly using the DVODE solver ( Brown 



et al. 1989). In cases where the chemical or cool- 



ing timesteps are much shorter than the hydrody- 
namic time step, sub-cycling is used to treat the 
cooling and chemistry, thereby avoiding the need 
to constrain the global time step. Chemical abun- 
dances are represented as tracer fields and are ad- 
vected using the standard FLASH infrastructure 
for such fields. 

2.3. Initial Conditions and Resolution 

Our initial conditions are motivated by typical 
ISM values. At t — 0, we start from a homoge- 
neous box with a total mass of 3.99 x 10^ AIq. The 
box length is L — 500 pc, the initial gas density 
po = 2.16 X 10~^^ g cm^"^, and the initial tem- 
perature is set to 5000 K. Turbulent driving with 
the energy input rate given in Section 2.1 results 



in an rms turbulent velocity of 50 km s""'^, and so 
the initial rms Mach number is approximately 7. 
We show results for two different gas metallicities: 
solar metalhcity (hereafter ZO) and low metalhcity 
Z = IQ-^Zq (hereafter Z-3). We use a 512^ grid 
without AMR for all of the runs presented. With 
this setup each cell has a volume of (0.977 pc)'^. 

3. Results 

We first evolve the box by driving turbulence 
for one turnover time Tdyn = 50 Myr. Thereafter, 
driving is switched off, and the gas is evolved for 
a further turnover time, up to t = 100 Myr. In 
the following sections, we discuss the most striking 
features and differences of these two evolutionary 
stages with respect to the gas metalhcity. 

3.1. Driven turbulence 

Starting from a homogeneous gas density distri- 
bution, a two-phase medium is quickly generated 
due to the infiuence of driven, supersonic turbu- 
lence, which produces strong non-linear fiuctua- 
tions in all thermodynamic variables. Kinetic en- 
ergy is injected on large scales and cascades to in- 
creasingly smaller ones. This turbulence generates 
shocks that compress the warm, neutral medium 
into dense, cold filaments. 

Regardless of the gas metalhcity, regions with 
very low density {p < 10^^''g cm^'^) and very high 
density (p > 4 x lO^^^g cm^"^) develop quickly, 
after only 5 Myr = 0.1 Tdyn- The metalhcity influ- 
ences the maximum density reached, Pmax, which 
is about a factor of eight larger in the solar metal- 
hcity run compared to the low-metallicity run. 
However, this is a relatively small difference com- 
pared to the total range of densities. 

In isothermal compressible turbulence the 
volume-weighted density PDF is expected to fol- 



low a log-normal distribution (Mac Low 



Ostriker et al . 2001; Kritsuk et al.|2007[|FedeTrath 



1999 



et al.||2008[lkitsionas'et al.||2009[ ) of the form 



p{s)ds 



: exp 



2a2 



ds. (1) 



where s = ln{p/pQ) is the logarithmic density con- 
trast, with pq equal to the mean gas density. In 
Equation ([T]), cr is the logarithmic density disper- 
sion, and sq is the mean value of s. 
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-5 5 1 3 3 4 

s = In (p/pa) s = logio(T) 



Fig. 1. — Density and temperature PDFs after 1 Tdyn of evolution under driven conditions (upper panels) 
and after another Tdyn of evolution under decaying conditions (lower panels). Solid black lines: case ZO; 
dashed black lines: case Z-3; green lines: log-normal fits (see Table 1). 
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In Figure [T] (upper left panel) we show the re- 
sulting density PDFs after driving turbulence for a 
turnover time (i.e., five turbulent crossing times). 
Where appropriate, we include fits according to 
Equation ([I]), and summarize the resulting values 
for cr and sq (or p{so) = /Ooexp(so), respectively) 
in Table 1. In the driven case, the density PDFs 
are well described with a single log-normal func- 
tion for both ZO and Z-3. For ZO the distribution 
is slightly asymmetric and skewed toward higher 
densities. However, the best fits result in a simi- 
lar variance of cr ~ 1.7 in both cases. Studies of 
driven turbulence in isothermal gas have found a 
relationship between a and the rms Mach number 



M ( 


Padoan et al.|1997 


Padoan & Nordlund|2002 


Price et al.||2011 


) 



ln(l 



(2) 



where the value of the parameter b depends on 
the relative importance of the compressive and 
solenoidal modes ( Federrath et al.pOOS ). For the 
driving used in our simulations, b = 0.5, and the 
rms Mach number at t = 50 Myr is 8.66 in the 
solar metallicity simulation, and 9.26 in the low- 
metallicity simulation. Equation ^ therefore pre- 
dicts that cr = 1.72 for the solar metallicity simu- 
lation and a ~ 1.76 for the low metallicity simula- 
tion, in good agreement with our measured value 
of 1.7. 

We emphasize that our finding of a broad, but 
single-peaked density PDF supports the conclu- 



sions of Vazquez- Semadeni (20091, where he ar 



gues that supersonic turbulence should dominate 
the evolutionary state of ISM gas because it acts 
faster than conventional TIs. We note that a 
single-peaked density PDF has also been found by 
de Avillez & Breitschwerdt (20071. However, mod- 



els by e.g., Slyz et al. (2005 1, who drive turbulence 



in colliding flows, result in a double-peaked den- 
sity PDF. We suggest that the transition between 
single- and double-peaked density PDF profile de- 
pends on the strength of the turbulence, i.e., the 
Mach number. This suggestion is verified below 
(Section 3.2). 

The temperature PDFs, which are also shown 
in Figure [l] (upper right panel) cannot be well ap- 
proximated with a single log-normal distribution, 
which only gives acceptable fits for the high tem- 
perature parts. Toward low temperatures (i.e.. 



for T < 300K), the distributions show extended 
powerlaw tails. This behavior is caused by the 
cooling instability that separates the warm and 
cold neutral gas. 

Overall, we find the temperature distributions 
to be more strongly affected by the gas metallicity 
than the density PDFs. In particular, the volume 
filling factors (VFFs) of the cold components differ 
significantly. For further analysis, we define the 
cold component to be all gas with T < BOOK, as 
300 K is the temperature threshold below which a 



stable, cold atomic phase can exist (Wolfire et al 



2003). Then we find 1.7% of the computational 
domain to be cold for ZO, whereas we find only 
0.075% to be cold for Z-3. 

In Figure [2j columns A and B, we show the 
density (top panels) and temperature (bottom 
panels) of the evolved filamentary structures in 
the z — plane. Prima facie the density distri- 
butions of the solar metallicity (panel (A)) and 
the low-metallicity case (panel (B)) look similar, 
in agreement with the PDFs. However, the fila- 
mentary structure formed in run Z-3 seems to be 
coarser and individual filaments seem to be less 
condensed than in run ZO. 

When comparing the two-dimensional density 
and temperature distributions (Figure [2| , we find 
both shock heating and pdV heating as well as 
radiative cooling and cooling by expansion to be 
very important. The presence of shock- heated re- 
gions surrounding dense and cool filaments can be 
found in both the ZO and the Z-3 case. However, 
most of the cold regions in run Z-3 can be identi- 
fied with voids in the density structures (see also 
Figure 3). Thus, these regions have cooled due 
to expansion. Radiatively cooled regions within 
dense filaments are very rare in run Z-3, whereas 
these dominate the cold part in run ZO. It should 
be noted that for ZO the mean density in the cold 
gas is higher than the mean density in the warm 
gas (T > 300K), as expected for gas which truly 
undergoes a cooling instability. However, for Z-3 
this is not the case (see Table 1). This indicates 
that cooling due to expansion is dominant in case 
Z-3 whereas the dense, cold component is very 
small. Thus, our results confirm previous semi- 



analytical estimates (e.g., Schaye 2004) that the 
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1 00 200 300 400 500 1 00 200 30O 400 500 100 200 300 400 500 100 200 300 400 500 
K [pc] It IpD] JC IM le IM 



Fig. 2. — Cuts through the computational volume on the (z = 0)-plane, showing the density and temperature 
of the developed two-phase media after supersonic turbulence has been driven for five turbulent crossing 
times (panels (A) and (B)), and then allowed to decay for a further five turbulent crossing times (panels 
(C) and (D)). Panel (A) and (C) depict the solar metallicity case, ZO, whereas panel (B) and (D) show the 
low-metallicity case, Z-3. 
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formation of stable molecular cloud cores might 
be strongly inhibited in low-metallicity gas. 



3.2. Decaying turbulence 

Starting from the evolved turbulent setup dis- 
cussed before, we switch off turbulent driving at 
Tdyn = 50 Myr, and study the evolution of the 
gas under decaying conditions until t = 2 Tdyn- 

In 

Figure 3 (upper panel) we show the evolution of 
the rms velocity dispersion, Vrms, for the full time 
interval. After driving has been switched off v^^is 
decreases. From initially close to 50kms^^, Vnns 
drops down to 20kms~^ within 5 Myr {t — 55 
Myr) and to lOkms^^ within approximately 25 
Myr {t ~ 75 Myr). In Figure 3 (lower panel) we 
depict the decay of kinetic energy. We find that 
both cases, ZO and Z-3, are well described by the 



fitting formula derived by Pavlovski et al. ( 2002 1 



E{t)=Eo{l + t/t,y\ (3) 

where Eq is the initial kinetic energy and ti = 
-^inj/vrms IS cqual to the initial ratio of mean 
energy injection scale to the rms velocity. For 
isothermal, compressible MHD turbulence, Mac 
Low (19991 derived rj w —1. Kitsionas et al. 



(2009) also extracted ry = — 1 for compressible HD 
turbulence. Incompressible, non-magnetic turbu- 
lence has been found to decay somewhat faster 
with T] in [—1.2,-1.4]. In addition, coohng has 
been found to accelerate the decay as the aver- 



age Mach number is larger (Kritsuk & Norman 



2002 Pavlovski et al. 2002 1. With cooling and 



MHD turbulence, [Pavlovski et al.| (|2002[) d erived 
T] « —1.3, whereas Kritsuk & Norman (2002) 



found indices between —1 and —2. We derive 
T] « —1.3 for Z-3. For ZO we find a slow decay 
within the first 10 Myr (77 w —1), but a steep de- 
cay thereafter with 77 w —1.8 from 10 Myr {t = 60 
Myr) onward. Due to the fact that cooling is much 
more efficient in the case of solar mctallicity, ZO 
(see Section 4.1 for a discussion of the mass frac- 
tion in cold gas), and based on previous findings 
that the decay is accelerated with cooling, our 
findings are in overall agreement with the expec- 
tation that r] should be lower for ZO than for Z-3. 

Returning to Figure 1, we show the tempera- 
ture and density structure of run ZO (panel (C)) 



and run Z-3 (panel (D)) after turbulence has de- 
cayed for 50 Myr. In this case we find striking 
differences between ZO and Z-3. Whereas cold 
clumps survive in run ZO, they are completely dis- 
persed in run Z-3. In fact, we find all of the cold 
gas in Z-3 to be warmer than 300 K after 10 Myr 
(i.e., at t = 60 Myr). 

In agreement with this qualitative finding, the 
metallicity also has a major impact on both den- 
sity and temperature PDFs (see lower panels in 
Figure [T]) . All distributions show a significant de- 
crease in variance as the system evolves toward 
its equilibrium density and temperature configu- 
ration. In run Z-3 a single-peaked density PDF 
remains. However, now the variance of the distri- 
bution is so small that we can speak of a single, 
warm phase. All of the cold gas has been heated 
and therefore the VFF of the cold component is 0. 
In this environment star formation is basically im- 
possible. This is in agreement with the findings of 



Schaye (2004) and Spaans & Norman (1997), who 



predict that star formation is strongly suppressed 
in low-metallicity environments. 

For run ZO, the single log-normal density PDF 
is replaced by a combination of two log-normals 
representing two distinct gas phases, a warm dif- 
fuse and a cold dense phase. A double-peaked 
density PDF has been previously reported by au- 



thors like Slyz et al. (2005) and Audit & Hen 



nebelle (20101, who employed cooling functions. 



In contrast to their work, the separation into two 
distinct phases seems to be less clean in our case. 
We note that the blurring of the warm and cold 
phases is probably a consequence of our chemo- 
dynamical model. We expect the separation to 
become increasingly clear as the turbulence de- 
cays further and Vrms decreases. We will discuss 
the time evolution of the density PDFs for dif- 
ferent turbulent driving energies in a subsequent 
paper. Finally, the VFF or of the cold component 
is slowly increasing in the case of ZO, and reaches 
4% at t = 100 Myr. 

4. Discussion 

In the following, we discuss the impact of tur- 
bulence and non-equilibrium chemistry on the dy- 
namical evolution of the ISM. 
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SS.M 72.50 50.00 S7,S0 5.00 86.00 75.50 50.00 S7,S0 5.00 




Enclosed mass [%] Enclosed mass [%] 



85,00 73.60 50.00 37.50 6.00 86.00 7S.50 50.00 S7,S0 5.00 




95,00 72.50 50,00 27,50 5,00 95,00 73.50 50.00 37.50 5.00 



Fig. 4. — Cumulative mass distributions derived for driven turbulence at t = 50 Myr (top row; see Figure 
2 columns A and B) and for decaying turbulence at t = 100 Myr (middle row; see Figure 2 columns C 
and D). The equilibrium curves of density and temperature for Zq = 1 Zq (black solid line, left column) 
and Zq = 10~^ Zq (black dashed line, right column) are plotted on top. For comparison we show the mass 
distributions of runs EI and KI (bottom row), both at the end of the driven period at t = 50 Myr. See 
Section 4.2 for further details on these runs i?/and KI. 
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Run 


Time 


Max 




Mean Warm 


Mean Cold 


Fit 






[Myr] 


[g cni~^ 


] 


[g cm-3] 


[g cm-3] 






<T 


ZO 


50 


2.4 X 10" 


21 


1.65 X 10^24 


3.35 X 10-23 


4.4 X 10- 


25 


1.73 


Z-3 


50 


3.2 X 10- 


22 


2.07 X 10-24 


7.33 X 10-25 


4.6 X 10- 


25 


1.70 


ZO 


100 


1.0 X 10~ 


21 


5.83 X 10^25 


1.22 X 10-23 


6.3 X 10- 


25 


0.22 














8.2 X 10- 


24 


0.93 


Z-3 


100 


1.2 X 10" 


23 


9.27 X 10-25 




1.9 X 10- 


24 


0.31 



Table 1: Absolute Maximum (Column 3) and Mean Densities (Columns 4 and 5) Extracted from the Simu- 
lations. The first two rows are for the part of driven turbulence, the following three rows are for the part of 
decaying turbulence. The mean densities are time averages over 50 Myr for the warm (T > 300K, Column 4) 
and cold (T < 300 K, Column 5) gas components individually. For reference, the mean gas density summed 
over both components was 2.16 x 10-2^ g cm-^ throughout. We also list the parameters for the log-normal 
fits to the density PDFs shown in Figure [l] The parameters are described in Equation ([T]) . For run ZO 
we list two sets of parameters &t t — 100 Myr as this distribution can only be fitted with the sum of two 
log-normal functions (see Figure 1). 



4.1. Metallicity dependence 

In Figure |4] (upper and central panels) we show 
the equilibrium density-temperature curves for 
low and solar metallicity together with the mass 
distributions found in ZO and Z-3 at 50 Myr and 
100 Myr. At low densities (p < lQ~'^^g cm-^) 
and high temperatures, cooling is dominated by 
Lya and therefore the metallicity dependence is 
negligible. For higher densities and lower tem- 
peratures, the cooling rate is dominated by fine 
structure emission from C^ and O, and is there- 
fore strongly dependent on Z. 
Based on comparing the equilibrium curves for 
low and high Z, it might at first seem surpris- 
ing that both runs {ZO and Z-3) look so simi- 
lar in case of driven turbulence and show a log- 
normal profile which is single-peaked, just like in 
an isothermal case. Therefore we conclude that 
the shape of the density PDF (single or double 
peaked; see Figure [T]) is not only a function of the 
employed gas thermodynamics but at the same 
time must be a function of the turbulent Mach 
number. Strongly driven turbulence (Vazquez 



Myr), where Wrms has decayed to 7 kms the 
separation into two phases is less distinct than in 



Semadeni|2009 ) causes the gas to populate a large 
region of parameter space in the (p, T)-plane sig- 
nificantly displaced from the equilibrium locus. 
As a result, the density PDF in the case of driven 
turbulence is broad, but single-peaked just like in 
an isothermal simulation (Figure [l] upper pan- 
els). Even at late stages of the evolution {t = 100 



previous models (Slyz et al. 2005 Audit & Hen- 
nebellepOlO ). A closer investigation reveals that, 
apart from the different thermodynamics schemes 
employed, also the Mach number seems to be lower 
in these models than in our simulations. For in- 
stance, Audit & Hennebelle (20101 who reported 



a clearly double peaked profile had generally low 
turbulent Mach numbers (the Mach number is 
around 1 in Audit & Hennebelle 2010). On the 



other hand, authors who used efficient driving 



mechanisms (de Avillez & Breitschwerdt 2007) 



also reported single-peaked distributions despite 
including gas heating and cooling. We therefore 
infer that a pronounced double-peaked density 
PDF is a sign of a two-phase medium at low Mach 
numbers. 

Even though the density PDFs in the driven 
turbulent cases look similar and seem to be in- 
dependent of the gas metallicity, we find the re- 
spective mass distributions in the (p, T)-plane to 
be quite different. For ZO, a large fraction of the 
gas is in the thermally unstable regime, in agree- 



ment with observations (e.g., Kanekar et al. 20031. 
Moreover, for ZO, low temperature values prefer- 
entially coincide with high-density values. For Z- 
3, the opposite trend is found. Here, lower tem- 
perature values correspond to lower densities. In 
general, only a small mass fraction is found to 
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50.5 



Fig. 3. — Upper panel: we depict the root mean 
square velocity as a function of time. From left to 
right, the first dotted vertical line marks the begin- 
ning of the decay phase at t = 50 Myr, the follow- 
ing two dotted horizontal and vertical lines mark 
the points at which the turbulence has decayed to 
Vrms = 20kms~'^ (55 Myr) and v^-^s = lOkms"^ 
(75 Myr), respectively. Lower panel: decay of 
kinetic energy for the second 50 Myr of evolu- 
tion. We show both cases, ZO (solid line) and Z-3 
(dashed line), as well as the decay curves derived 
from Equation (3) (dotted lines) for -q = —1.3 and 
r] = —1.8. The kinetic energy decay is faster in 
solar metallicity gas. 




20 40 60 80 100 
Time [Myr] 



Fig. 5. — Time evolution of the volume filling 
factor (solid lines) and the mass fraction (dashed 
lines) of the cold gas (T < 300 K). Turbulent driv- 
ing was switched off at t = 50 Myr. We show all 
four runs: ZO (black), Z-3 (green), EI (turquoise), 
and KI (magenta). 



be at temperatures less than 1000 K, as cold gas 
preferentially coincides with regions of low density, 
where expansion is the main cooling process. 

We now determine the fraction and the VFF 
of the cold gas formed in our simulations. For 
decreasing metallicity the equilibrium tempera- 
ture, Toqui, remains high (Tcqui > 300 K) up 
to increasingly high density. For further refer- 
ence, we define the threshold density pthrcs above 
which gas may cool below 300 K. For Zo = Zq, 
we find pthros — 9 x lO"^** g cm~'^, and for 
Zq = 10~^Zq we find pthrcs ^ 5 x 10^^^ g cm~^. 
Thus, for Zq = Zq some of the gas reaches densi- 
ties higher than pthres, whereas for Zq = 10"^^©, 

Pmax(^-3) = 3.2xl0"2^g Cm-3 andpthros > Pmax 

(see Section 3.1 and Table 1). 



The strong influence of the gas metallicity on 
the formation and survival of cold coherent clumps 
only reveals itself during the non-driven, decaying 
phase, as the temperature and density PDFs of the 
two simulations are very similar during the driven 
turbulent period. After turbulent driving has been 
switched off (Figure |4j middle panels) the systems 
evolve toward local equilibrium. Consequently, 
the almost immediate disappearance of the cold 



11 



component in run Z-3 can be explained. The cold 
component formed in the driven case was solely 
a product of turbulent stirring. Regions of high 
density, cold gas, which truly undergo thermal in- 
stability are rare at low metallicity since the high 
threshold density required for the formation of a 
stable cold component cannot be reached. There- 
fore, when evolved under decaying conditions, the 
cold gas is quickly heated up and no cold and dense 
clumps remain. 

This result is reinforced by the time evolution 
of the VFF of the cold gas component, and of the 
mass in cold gas, shown in Figure [5j With driven 
turbulence the VFF is roughly constant and about 
20 times lower in Z-S than in ZO. With decaying 
turbulence the two cases diverge as VFF is slowly 
increasing for ZO whereas it quickly approaches 
for Z-3. The mass fraction of gas found in the 
cold component evolves in a similar manner as the 
VFF. However, even in the driven case, the mass 
in cold gas is more than a factor of 1000 lower 
for Z-3 than for ZO (with an average of 10^ Mq 
for Z-3, corresponding to only 0.025% of the total 
mass, compared to an average of 10^ Mq for ZO, 
corresponding to a substantial fraction of the total 
mass, i.e., 25%). 

4.2. The impact of non-equilibrium chem- 
istry 

In the following, we analyze the impact of non- 
equilibrium chemistry on our findings. To quantify 
the degree of any departures from thermal equi- 
librium in runs ZO and Z-3 we first monitor the 
mean difference between actual gas temperature 
and expected equilibrium temperature in Figure [6] 
(solid and dashed line). In the case of solar metal- 
licity, ZO, we find a significant offset, AT, of the 
actual gas temperature from the equilibrium tem- 
perature, in agreement with Figure 4 (upper left 
panel). In the decaying period (t > 50Myr) AT 
then decreases, as does the turbulent Mach num- 
ber. For Z-3, the system is close to equilibrium 
throughout the simulation and AT is small. 

Two possible causes of the high AT in run ZO 
are: (1) non-equilibrium chemical effects and (2) 
efficient thermalization of turbulent energy on a 
time scale shorter than the cooling time scale. To 
investigate the importance of the non-equilibrium 
chemistry on the gas heating and cooling, we per- 
formed two comparison runs at Z = Zq that do 



not take into account non-equilibrium effects. In 
run EI, we still consider the detailed Glover et al. 
cooling function when solving the energy equation, 
but at the same time we assume chemical equilib- 
rium. In run KI we completely neglect the chem- 
ical network as well as the detailed cooling func- 
tion, and instead use the parameterized cooling 
function: 



F 

A(T) 

r 



2.0 X 10-2^er" -"^ 



= 10^ cxp 



-\- 1.4x10" 



■g s 
-1.184 X 10' 
T + 1000 



Tcxp 



-92 



(4) 



(5) 



which is based on Koyama & Inutsuka ( 2000 2002 ) 
(see also Vazquez-Semadeni et al.||2007 |. 

We include the results for both EI and KI in 
Figures |4] (bottom row) , [5] and [6] From Figure 
|4j we find that the phase space mass distribution 
formed in run EI is almost identical to run ZO. 
Thus, non-equilibrium chemistry seems to have a 
minor effect at high driving velocities. We con- 
clude that the offset of the mass distribution with 
respect to the equilibrium curve as well as the high 
AT in Figure[6]are caused by turbulent shock com- 
pression and extra heating associated with turbu- 
lent dissipation. Run KI supports this picture, 
although we find KI to differ from EI and ZO in 
mass distribution (see Fig. [4j bottom right panel) 
and mean temperature (see Figure [6]). Overall, 
run KI is closer to thermodynamical equilibrium 
throughout the whole simulated time. 

In order to elaborate the global effect of using a 
detailed cooling function together with a chemical 
network rather than a simple, parameterized cool- 
ing function, we calculated the time evolution of 
the cold mass fractions and the VFFs for EI and 
KI and included them in Figure [5] (turquoise and 
magenta lines). We find that run EI is very sim- 
ilar to ZO throughout the driven phase, although 
cooling seems to be slightly stronger in case of ZO, 
as the VFF of the cold component is somewhat 
higher with respect to EI. On the other hand, in 
the decaying phase, this effect is reversed after 
t = 70 Myr, where frms ~ 13kms~^. The reduced 
cooling occurring in run ZO with full chemistry is 
also visible in the cold gas fraction, which drops 
below EI after t = 70 Myr. 
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In summary, we cannot see strong immediate 
effects of the non-equilibrium chemistry. We note 
that this finding might have to be revised once 
the formation of molecules (e.g., H2) is consid- 
ered (Gnedin & Kravtsov 2011). Furthermore, 



our simulations show that using a detailed, time- 
dependent cooling prescription is still superior 
to using a parameterized cooling function as the 
results are much more accurate. Run KI continu- 
ously overpredicts the fraction of cold gas (by up 
to 200%) and the VFF of the cold gas. 




0.0 



so 



40 60 
Time [Myr] 



SO 



mo 



Fig. 6. — Mean deviation of actual gas temper- 
ature and equilibrium temperature as a function 
of time, again for ZO (solid line) and Z-3 (dashed 
line). In addition, we show run EI (dash-dotted 
line) and run KI (dash-dot-dot-dotted line) . 



4.3. Other effects that may influence the 
phase distribution 

As the variance of the density PDF has been 
shown to be directly related to the turbulent Mach 



number ( ] 


^adoan et al. 


1997 


Padoan & Nord- 


lund 


2002 


Federrath et al. 


2008 1, the spread in 



p about the given equilibrium configuration is ex- 
pected to increase with increasing Mach number 
M. More compressive driving would enhance this 
effect ( Federrath et al.|["2010 |. Very strong turbu- 
lent driving could therefore result in the formation 
of a stable cold phase even if the metallicity were 
low. Also our models do not include the effects 
of self-gravity. Self-gravity is known to cause a 
power-law tail toward high densities in the den- 



sity PDF (Scalo et al. 


1998 


Klessen et al. 


2000 


Wada & Norman 2001 


Li et al. 2003| 


Slyz et al.| 


2005 Kritsuk et al. 2011 Cho & Kim 


2011) and 



therefore would certainly influence the gas density 
distribution at high p. As the powerlaw tail is 
typically flatter than the log-normal distribution 
(with power spectral indices of around —1.5) more 
dense gas would be available for cooling, and the 
cold component would become more massive lead- 
ing to more star formation. Therefore, the effect 
of self-gravity might shift the critical values of tur- 
bulent Mach number above which the formation 
of cold gas is suppressed, and the critical values of 
metallicity below which the formation of cold gas 
is suppressed. 

The assumed FUV heating rate, Fpuv, also in- 
fluences the distribution of phases. In our present 
study, we assume that the photoelectric heating 
rate per unit dust mass is the same in both simu- 
lations, so that the only change in Fpuv between 
the two simulations comes from the change in the 
total dust mass, which we assume to scale linearly 
with metallicity. In reality, the composition and 
size distribution of dust grains in very low metal- 
licity gas may be quite different from what we are 
used to in the Milky Way, which would lead to an 
additional change in Fpuv- The strength of the 
FUV radiation field is also an influential uncer- 
tainty. Decreasing the FUV intensity by a factor 
of 100 results in a decrease in pthres of about a 
factor of 10 ( Wolfire et al.|[2003 ), thus facilitating 
the formation of a cold component in regions of 
low star formation. In an upcoming paper we will 
further discuss the relationship between pthrcsi 
and FpuVi as well as that between VFF and M. 

Placed in an astrophysical context, our finding 
implies that star formation would effectively be 
inhibited in a system of low metallicity if no event 
capable of driving strong turbulence (e.g., spiral 
shocks in a disk galaxy, a merger, or the infall of 
gas onto a galactic disk) occurred. 

5. Conclusions and Astrophysical Context 

The dependence of the cooling rate on the metal 
abundance in the gas shifts the pressure range that 
allows for the formation of a bi-stable medium 
toward higher values with decreasing metallicity. 
Even in the Milky Way, this effect is so influential 
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that the formation of cold molecular gas is pre- 
vented at distances larger than 18 kpc from the 
Galactic Center where Z k, 0.2 Zq (see analyti- 
cal estimates by |Wolfire et al. 12003) |Elmegreen fc| 
Elmegreen l[T984l ). 

In order to understand this process in greater 
detail, we have simulated the turbulent evolution 
of a multi-phase ISM under driven and decay- 
ing turbulent conditions at different metallicities 
(Z = Zq and Z = IO^^Zq). This mimics in- 
terstellar conditions for different galactic environ- 
ments. For driven turbulent conditions, we find 
that metallicity does not have a strong influence 
on either the shape of the resulting, single-peaked 
density PDF, or the temperature PDF. However, 
for Z = 10^"^ ^0 the VFF of the cold component is 
lower by an approximately constant factor of 20. 

In decaying turbulent conditions, we find a 
clear relation between the abundance of cold gas 
and the metallicity. Whereas the coldest regions 
are quickly dispersed in case of low metallicity 
and the system evolves toward a single compo- 
nent, warm medium, a stable cold component is 
formed in case of solar metallicity and a double- 
peaked density PDF is developed. Thus, on a 
dynamical time scale, the survival of cold clumps 
is strongly influenced by the gas chemistry. In 
agreement with previous models of other authors 
(e.g., Schaye|2004 1 , our numerical simulations im- 
ply that the turbulent evolution of the ISM at 
different metallicities causes systematic changes 
in the initial conditions for star formation and 
should strongly affect the global star formation 
efficiency of a galaxy. 

We also find the non-equilibrium chemistry to 
not have a significant influence on our results at 
high turbulent velocity dispersions for both, low 
and solar metallicity. However, in mildly turbulent 
conditions (wrms ^ 13 km s~^), the effect of non- 
equilibrium chemistry becomes more pronounced. 
Moreover, we find that using an accurate cool- 
ing prescription, like the one presented in |Glover| 
& Mac Low (2007), rather than a parameterized 



cooling function is especially important to pre- 
cisely predict the mass and VFF of the cold gas 
component in a multi-phase model. 

Placed in a more global astrophysical context 



these results could be meaningful to several fields. 
We can speculate on the following. 

• Metal-poor dwarf galaxies have very low star 
formation rates ( Boissier et al.|[2003 ), which 
can currently not be reproduced in simula- 
tions of dwarf galaxy formation and evolu- 



tion 


Quo et al. 


2011 


Gnedin & Kravtsov 


2010 


. If an external trigger (e.g., a merger 



or tidal perturbation) were missing, our 
chemo-dynamical ISM model would lead to 
such a low star formation rate. 

• Dwarf galaxies seem to have complex star 
formation histories (e.g., as observed for the 
Carina dSph by |Tolstoy| [20031 [Burkert fc| 
Ruiz-Lapuente 19971. In our model, metal- 



poor dwarfs could only form stars if the 
gas were stirred. Hence, clumpy gas infall, 
tidal effects, or other forms of irregular stir- 
ring (e.g., due to gravitational instabilities) 
would naturally lead to very complex star 
formation histories in these systems. 

• Metal-poor galaxies which have a low star 
formation rate should be able to slowly ac- 
crete a lot of gas - e.g., via gas infall in cold 
streams ( [Keres et al.|2005||Dekel et al.|2009] ) 
- without experiencing a significant increase 
in their star formation rate. Therefore, they 
would acquire the potential to become pro- 
genitors for massive starburst galaxies. 

• The star formation rate in merging galaxies 
is probably not strongly dependent on the 
gas metallicity if supersonic gas turbulence 
is driven due to the merger. 
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